// -*- C++ -*-----------------------------------------------------------------
// 
// Project name:    BiocompUP Align
//
// Date:            02/2001
//
// 
// -----------------x-------------------x-------------------x-----------------

#include <string>
#include <GetArg.h>
#include <AlignmentBase.h>
#include <StructuralAlignment.h>
#include <PdbLoader.h>
#include <PdbSaver.h>

using namespace Biopool;

void sShowHelp()
{
  cout << "CE to (structural) Alignment \n"
       << "\n"
       << " Options: \n"
       << "\t-i <filename> \t\t Input CE file\n"
       << "\t--pdb1 <filename> \t\t First PDB file (''target'')\n"
       << "\t--pdb2 <filename> \t\t Second PDB file (''template'')\n"
       << "\t[-m <filename>] \t\t Max CA distance for ''equivalence''\n"
       << "\t[-o <filename>] \t\t Output PDB file\n"
       << "\n";
};

void sAddLine()
{
  cout << "-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-\n";
}


int main(int nArgs, char* argv[])
{ 
  // --------------------------------------------------
  // 0. treat options
  // --------------------------------------------------

  if (getArg( "h", nArgs, argv))
    {
      sShowHelp();
      return 1;
    };

  string inputFile, outputFile, pdbFile1, pdbFile2;
  double maxDist;
  getArg( "i", inputFile, nArgs, argv, "!");
  getArg( "-pdb1", pdbFile1, nArgs, argv, "!");
  getArg( "-pdb2", pdbFile2, nArgs, argv, "!");
  getArg( "o", outputFile, nArgs, argv, "!");
  getArg( "m", maxDist, nArgs, argv, 4.0);

  if ((inputFile == "!") || (pdbFile1 == "!") || (pdbFile2 == "!"))
    {
      cout << "Missing file specification. Aborting. (-h for help)" << endl;
      return -1;
    }

  // --------------------------------------------------
  // 1. read alignment
  // --------------------------------------------------

  StructuralAlignment ali;

  ifstream aliFile(inputFile.c_str());
  if (!aliFile)
    ERROR("File not found.", exception);

  ifstream pdbIn1(pdbFile1.c_str());
  if (!pdbIn1)
    ERROR("File not found.", exception);
  PdbLoader pl1(pdbIn1);

  Spacer spTarget;
  spTarget.load(pl1);

  ifstream pdbIn2(pdbFile2.c_str());
  if (!pdbIn2)
    ERROR("File not found.", exception);
  PdbLoader pl2(pdbIn2);
  pl2.setNoHAtoms();
  pl2.setNoConnection();

  Spacer spTemplate;
  spTemplate.load(pl2);

  ali.setTarget(spTarget);
  ali.loadCE(aliFile, spTemplate);

 

  sAddLine();

  ali.buildEquivalenceNetwork();
  ali.buildFragmentNetwork(maxDist);
  
  ali.writeData(maxDist);
  
  sAddLine();

  cout << "B-factor comparison for fragments:\n";

  for (unsigned int i = 0; i < ali.fragData.size(); i++)
  {
      double Bfac1 = 0.0;
      for (unsigned int j = ali.fragData[i].min; j <= ali.fragData[i].max; j++)
	  if (j < ali.getTarget().sizeAmino())
	      Bfac1 += ali.getTarget().getAmino(j)[CA].getBFac();
        Bfac1 /= ali.fragData[i].max - ali.fragData[i].min + 1;

      double Bfac2 = 0.0;
      for (unsigned int j = ali.fragData[i].min; j <= ali.fragData[i].max; j++)
	  if (j < ali.getTemplate().sizeAmino())
	      Bfac2 += ali.getTemplate().getAmino(j)[CA].getBFac();
        Bfac2 /= ali.fragData[i].max - ali.fragData[i].min + 1;

      cout << "# " << i << "\t" << setw(5) << setprecision(3) << Bfac1
	   << "\t" << setw(5) << setprecision(3) << Bfac2 << "\n";

      // WARNING: indexes for B-factor calculation appear to be shifted
      // compared to the original PDB file!!
  }

  sAddLine();

  // --------------------------------------------------
  // 3. write alignment to disk
  // --------------------------------------------------

  if (outputFile != "!")
  {
      ofstream outFile(outputFile.c_str());
      if (!outFile)
	  ERROR("File not found.", exception);
      PdbSaver ps(outFile);

      ali.getTarget().save(ps);
      ali.getTemplate().save(ps);
      ps.endFile();
  }
      
}
